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Two-photon processes have so far been considered only as resulting from frequency-matched second-order 
expansions of light-matter interaction, with consequently small coupling strengths. However, a variety of novel 
physical phenomena arises when such coupling values become comparable with the system characteristic fre¬ 
quencies. Here, we propose a realistic implementation of two-photon quantum Rabi and Dicke models in 
trapped-ion technologies. In this case, effective two-phonon processes can be explored in all relevant parameter 
regimes. In particular, we show that an ion chain under hichromatic laser drivings exhibits a rich dynamics and 
highly counterintuitive spectral features, such as interaction-induced spectral collapse. 


I. INTRODUCTION 

The quantum Rabi model describes the interaction of a two- 
level quantum system, a qubit, with a quantized single-mode 
bosonic field. Its semiclassical limit, where a classical field 
is considered, is known as the Rabi model [1]. In the last 
few decades, the quantum Rabi model has been used in a 
regime where the rotating-wave approximation (RWA) holds, 
giving rise to the Jaynes-Cummings model [2] and describing 
a plethora of experiments, mostly related to cavity quantum 
electrodynamics. On the other hand, from a mathematical 
point of view, an analytical solution for the spectrum of the 
quantum Rabi model has been recently developed [3]. Such 
results have prompted a number of theoretical efforts aimed 
at applying similar techniques to generalizations of the quan¬ 
tum Rabi model, including anisotropic couplings [4-6], two- 
photon interactions [7, 10] and multi-qubit extensions [11], as 
is the case of the Dicke model. 

In particular, the two-photon quantum Rabi model enjoys a 
spectrum with highly counterintuitive features [12, 13], which 
appear when the coupling strength becomes comparable with 
the bosonic mode frequency. In this sense, it is instructive to 
compare these features with the ultrastrong [14—16] and deep 
strong [17] coupling regimes of the quantum Rabi model [18]. 
The two-photon Rabi model has been applied as an effective 
model to describe second-order processes in different physical 
setups, such as Rydberg atoms in microwave superconducting 
cavities [19] and quantum dots [20, 21]. However, the small 
second-order coupling strengths restrict the observation of a 
richer dynamics. 

In trapped-ion systems [22, 23], it is possible to control the 
coherent interaction between the vibrations of an ion crystal 
and internal electronic states, which form effective spin de¬ 
grees of freedom. This quantum technology has emerged as 
one of the most promising platforms for the implementation 
of quantum spin models, including few [24] or hundreds [25] 
of ions. Disparate complex quantum phenomena have been 
explored using trapped-ion setups, such as Ising spin frustra¬ 
tion [26], quantum phase transitions [27, 28] and the inhomo¬ 


geneous Kibble-Zurek mechanism [29]. Eurthermore, second 
sidebands have been considered for laser cooling [30] and for 
generating nonclassical motional states [31-33]. 

In this paper, we design a trapped-ion scheme in which the 
two-photon Rabi and two-photon Dicke models can be realis¬ 
tically implemented in all relevant regimes. We theoretically 
show that the dynamics of the proposed system is character¬ 
ized by harmonic two-phonon oscillations or by spontaneous 
generation of excitations, depending on the effective coupling 
parameter. In particular, we consider cases where complete 
spectral collapse—namely, the fusion of discrete energy lev¬ 
els into a continuous band—can be observed. 


II. THE MODEL 

We consider a chain of N qubits interacting with a single 
bosonic mode via two-photon interactions 

n n 

where h = 1, a and a) are bosonic ladder operators; cr" and 
cr" are qubit Pauli operators; parameters w, w”, and gn, repre¬ 
sent the mode frequency, the n-th qubit energy spacing and the 
relative coupling strength, respectively. We will explain below 
how to implement this model using current trapped-ion tech¬ 
nology, considering in detail the case = 1 and discussing 
the scalability issues for > 1. 

We consider a setup where the qubit energy spacing, Wint, 
represents an optical or hyperfine/Zeeman internal transition 
in a single trapped ion. The vibrational motion of the ion is de¬ 
scribed by bosonic modes a, with trap frequency iv. Turn¬ 
ing on a hichromatic driving, with frequencies and Wb, an 
effective coupling between the internal and motional degrees 
of freedom is activated. In the interaction picture, the standard 
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Hamiltonian [22] describing this model reads 

^ ^ + H.c.j , 

j=r,b 

( 2 ) 

where a{t) = a Here, 12,. and Hb are coupling param¬ 

eters directly proportional to the driving laser amplitude, and 
(pj represents the phase of each laser with respect to the atom 

dipole. The Lamb-Dicke parameter rjj = ^ defined 

by the the projection of the j-th laser field wavevector in 
the z direction and by the ion mass m. We consider the sys¬ 
tem to be in the Lamb-Dicke regime, 77 ^( 2 (fi)+l) ^ 1, where 
fi = a) a is the phonon number operator. 

We set the frequencies of the bichromatic driving to be 
detuned from the second sidebands, ojr = Wint — Iv + 6r, 
ujb = Wint + 2v + Sb- We choose homogeneous Lamb-Dicke 
parameters rjj = 77 , phases 4>j — 0 , and coupling strengths 
= 12 for both sideband excitations. Expanding the ex¬ 
ponential operator in Eq. (2) to the second order in p, and 
performing a RWA with 6j , 12^ ^ ly, we can rewrite the inter¬ 
action picture Hamiltonian 


n' = - 


77^12 ■ 


a2 ^-iS,t 


<7+ + H.c. (3) 


The first-order correction to approximations made in deriving 
Eq. (3) is given by j e^*^‘^*(T_|_-|-H.c., which produce spurious 

excitations with negligible probability Pg = (|^)^. Eurther 
corrections are proportional to pH or and oscillate at fre- 

quency ly, yielding Pe = ■ Hence, they are negligible 

in standard trapped-ion implementations. The explicit time 
dependence in Eq. (3) can be removed by going to another in¬ 
teraction picture with "Ho = | {Sb — 6r) a^a -b j (i^h -f Sr) (Tz, 
which we dub the simulation picture. Then, the system Hamil¬ 
tonian resembles the two-phonon quantum Rabi Hamiltonian 


o / i" f ‘2 ^ \ 

U^« = uj a'a+-Uz-g (Jx\a +a' J, 


(4) 


where the effective model parameters are linked to physical 
variables through u = \{Sr — Sb), Wq = —^{Sr + Sb), and 

g = Remarkably, by tuning Sr and Sb, the two-phonon 
quantum Rabi model of Eq. (4) can be implemented in all 
regimes. Moreover, the A-qubit two-phonon Dicke model of 
Eq. (1) can be implemented using a chain of N ions by apply¬ 
ing a similar method. In this case, the single bosonic mode 
is represented by a collective motional mode [37] (see Ap¬ 
pendix A). 

The validity of the approximations made in deriving Eq. (4) 
has been checked comparing the simulated two-photon quan¬ 
tum Rabi dynamics with numerical evaluation of the simulat¬ 
ing trapped-ion model of Eq. (2), as shown in Eig. 1. Stan¬ 
dard parameters and dissipation channels of current setups 
have been considered. In all plots of Eig. 1, the vibrational 
frequency is v j2-7 = 1 MHz and the coupling coefficient is 
H/27r = 100 KHz. The Lamb-Dicke parameter is p = 0.04 
for Eig. la and Eig. lb, while p = 0.02 for Eig. Ic. Notice 
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FIG. 1. (color online) Real-time dynamics for N = 1, resonant qubit 
LOq = 2uj, and effective couplings: (a) g = O.Olw, (b) g = 0.2u), 
and (c) g = 0.4a;. The initial state is given by \g, 2), i.e., the two- 
phonon Fock state and the qubit ground state. In all plots, the red 
solid line corresponds to numerical simulation of the exact Hamilto¬ 
nian of Eq. (1), while the blue dashed line is obtained simulating the 
full model of Eq. (2), including qubit decay — Is, pure dephas¬ 
ing t 2 = 30ms and vibrational heating of one phonon per second. 
In each plot, the lower abscissa shows the time in units of uj, while 
the upper one shows the evolution time of a realistic trapped-ion im¬ 
plementation. In panel (c), the full model simulation could not be 
performed for a longer time due to the fast growth of the Hilbert- 
space. 


that larger coupling strengths imply a more favourable ratio 
between dynamics and dissipation rates. Hence, the imple¬ 
mentation accuracy improves for large values of g/ui, which 
correspond to the most interesting coupling regimes. 
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III. REAL-TIME DYNAMICS 

Depending on the ratio between the normalized coupling 
strength g and the mode frequency oj, the model of Eq. (1) 
exhibits qualitatively different behaviors. Two parameter 
regimes can be identified accordingly. For the sake of simplic¬ 
ity, we will consider the homogeneous coupling case = g, 
uj'^ = ujq, for every n, and we will focus in the resonant or 
near-resonant case ojq « 2 w. 

In accordance with the quantum Rabi model, we define the 
strong coupling (SC) regime by the condition g/oj ^ 1. Un¬ 
der this restriction, the RWA can be applied to the coupling 

terms, replacing each direct interaction g J2n 
with g’ where we defined the rais¬ 
ing/lowering single-qubit operators crj = (cr” ± ) / 2 . 

When the RWA is valid, the system satisfies a continuous 
symmetry, identified by the operator C = a^a -f 2 cr"cr", 
which makes the model superintegrable [3]. In the SC regime, 
the interaction leads to two-photon excitation transfers be¬ 
tween the bosonic field and the qubits, as shown in Fig. la. 
Jaynes-Cummings-like collapses and revivals of population 
inversion are also expected to appear [38, 39]. 

As the ratio g/uj increases, the intuitive dynamics of the 
SC regime disappears and excitations are not conserved (see 
Fig. lb and Fig. Ic). When the normalized coupling ap¬ 
proaches the value g ~ Q.loj, the RWA cannot be performed, 
and the full quantum Rabi model must be taken into account. 
We define the ultrastrong coupling (USC) regime as the pa¬ 
rameter region for which 0.1 < g/uj < 0.5. An analytical so¬ 
lution for the system eigenstates has been derived in [7]. How¬ 
ever, this approach relies basically on a numerical instability 
related to the presence of dominant and minimal solutions of 
an associated three-term recurrence relation [ 8 ] and gives no 
qualitative insight into the behavior of the spectrum close to 
the collapse point. While continued-fraction techniques are 
applicable in principle [ 8 ], only a few low-lying levels can 
be computed and the method fails again in approaching the 
critical coupling (see below). While the G-function derived 
in [9] allows for the desired understanding of the qualitative 
features of the collapse, its mathematical justification is still 
incomplete. On the other hand, direct numerical simulation 
becomes challenging close to collapse due to the large num¬ 
ber of excitations involved. Especially the dynamics of the 
two-photon Dicke model is demanding for classical numeri¬ 
cal techniques. 

In the SC/USC transition, the continuous symmetry ( 
breaks down to a Z 4 discrete symmetry identified by the op¬ 
erator 

N 

n = (- 1 )^ cr” exp ■ (5) 

n—1 

We will call 11 the generalized parity operator, in analogy 
with the standard quantum Rabi model [3]. Four invari¬ 
ant Hilbert subspaces are identified by the four eigenvalues 
A = {1,-1,i,—i} of n. Hence, for any coupling strength, 
the symmetry H restricts the dynamics to generalized-parity 


' N=1I n = 1 : 

|g.o>- 
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FIG. 2. (color online) (a) Generalized-parity chains for N = 1, 2. 
For simplicity, for A = 2, only one generalized-parity subspace is 
shown, (b) Quantum state fidelity between the system state \(l>{t)) 
and the target eigenstates \'tpn} during adiabatic evolution. The 
Hamiltonian at f = 0 is given by Eq. (1) with N = 1, cOq/oj = 1.9 
and g = 0. During the adiabatic process, the coupling strength is 
linearly increased until reaching the value g/u) = 0.49. For the blue 
circles, the initial state is given by the ground state — 0)) = 
I For black crosses, |(/)(f = 0)) = while for the red 

solid line, = 0)) = The color code indicates gener¬ 

alized parity as in Fig. 3a. Notice that, due to generalized parity 
conservation, the fourth excited eigenstate of the decoupled 

Hamiltonian is transformed into the third one I'i/ig) of the full Hamil¬ 
tonian. 

chains, shown in Fig. 2a, for N = 1,2. 

When the normalized coupling g approaches g = ujj2 (see 
Fig 3c), the dynamics is dominated by the interaction term and 
it is characterized by photon production. Finally, when g > 
10/2, the Hamiltonian is not bounded from below. However, 
it still provides a well defined dynamics when applied for a 
limited time, like usual displacement or squeezing operators. 


IV. THE SPECTRUM 

The eigenspectrum of the Hamiltonian in Eq. (1) is shown 
in Figs. 3a and 3c for N = 1 and N = 3, respectively. Dif¬ 
ferent markers are used to identify the generalized parity H of 
each Hamiltonian eigenvector, see Eq. (5). In the SC regime, 
the spectrum is characterized by the linear dependence of the 
energy splittings, observed for small values of g. On the con¬ 
trary, in the USC regime the spectrum is characterized by level 
crossings known as Juddian points, allowing for closed-form 
isolated solutions [13] in the single-qubit case. 

The most interesting spectral features appear when the nor¬ 
malized coupling g approaches the value u>/2. In this case. 
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the energy spacing between the system eigenenergies asymp¬ 
totically vanishes and the average photon number for the first 
excited eigenstates diverges (see Fig. 3b). When g = w/2, 
the discrete spectrum collapses into a continuous band, and 
its eigenfunctions are not normalizable (see Appendix B). 
Beyond that value, the Hamiltonian is unbounded from be¬ 
low [12, 13]. This can be shown by rewriting the bosonic 
components of Hamiltonian of Eq. (1) in terms of the effec¬ 
tive position and momentum operators of a particle of mass 

TO, defined as a; = {a + a'^) andp = {a — a'^). 

Therefore, we obtain 

r ^2 ~\ 

77 ? LU ^ T)^ ^ 

n = — (w - 2g S^) +{uj + 2g S^)x'^ 

2 I m^cu^ 

n 

where Sx ^ ^x- Notice that Sx can take values in¬ 

cluded in the interval {Sx) G [—1,1]. Hence, the parameter 
(w -f 2g) establishes the shape of the effective potential. For 
g < a;/2, the particle experiences an always positive quadratic 
potential. For g = uj/2, there are qubit states which turn the 
potential flat and the spectrum collapses, like for a free par¬ 
ticle (see Appendix B). Finally, when p > a;/2, the effective 
quadratic potential can be positive, for {Sx) < —a;/2p, or 
negative, for {Sx) > oj/2g. Therefore, the Hamiltonian (6) 
has neither an upper nor a lower bound. 


V. MEASUREMENT TECHNIQUE 

A key experimental signature of the spectral collapse (see 
Fig. 3a) can be obtained by measuring the system eigenener¬ 
gies [40] when g approaches 0.5a;. Such a measurement could 
be done via the quantum phase estimation algorithm [41]. A 
more straightforward method consists of directly generating 
the system eigenstates [42] by means of the adiabatic proto¬ 
col shown in Fig. 2b. When g = 0, the eigenstates 
of Hamiltonian in Eq. (1) have an analytical form and can 
be easily generated [43]. Then, adiabatically increasing g, 
the eigenstates \ijj^) of the full model can be produced. No¬ 
tice that generalized-parity conservation protects the adiabatic 
switching at level crossings (see Fig. 2b). 

Once a given eigenstate has been prepared, its energy can 
be inferred by measuring the expected value of the Hamil¬ 
tonian in Eq. (1). We consider separately the measurement 
of each Hamiltonian term. The measurement of cr" is stan¬ 
dard in trapped-ion setups and is done with fluorescence 
techniques [22]. The measurement of the phonon number 
expectation value was already proposed in Ref. [44]. No¬ 
tice that operators tr" and a commute with all transforma¬ 
tions performed in the derivation of the model. The expec¬ 
tation value of the interaction term ga!^ can be 

mapped into the value of the first time derivative of (cr") at 
measurement time t = 0, with the system evolving under 

T-Lm = oja^a -f ^cr” — ga^ -f This Hamiltonian 




g/(o 

FIG. 3. (color online) Spectral properties of the Hamiltonian Eq. (1), 
in units of uj, for tOq = 1.9, as a function of the coupling strength g. 
For g > 0.5, the spectrum is unbounded from below, (a) Spectrum 
for N = 1. Different markers identify the generalized parity of each 
eigenstate: green circles forp = 1, red crosses forp = i, blue stars 
forp = — 1, and black dots forp = —i. (b) Average photon number 
for the ground and first two excited states, for N — 1. (c) Spectrum 
for N — 3. For clarity, the generalized parity of the eigenstates is 
not shown. 

is composed of a part A = a;o^a -F ^cr" which commutes 
with ct”, [A, a!)] = 0, and a part B = -F 

which anti-commutes with ct", {B, cr"} = 0, yielding 

„ ^^i{A+B)t^-i{A-B)t^n^ -pfie 

time derivative of this expression at f = 0 is given by 
{[i{A +B) — i{A—B)]a'^) = 2i(Bcr"), which is proportional 
to the expectation value of the interaction term of Hamiltonian 

inEq. (4), = 2{ga^ (a^ + a^^^). 

The evolution under Hamiltonian T-Lm in the simulation pic¬ 
ture is implemented in the same way as the Hamiltonian in 
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Eq. (4), but selecting the laser phases (pj to be |. Moreover, 
expectation values for the generalized-parity operator If of 
Eq. (5) can be extracted following the techniques described 
in Appendix C. 


VI. DISCUSSION 

We have introduced a trapped-ion scheme which allows one 
to experimentally investigate two-photon interactions in un¬ 
explored regimes of light-matter coupling, replacing photons 
in the model by trapped-ion phonons. It provides a feasible 
method to observe an interaction-induced spectral collapse in 
a two-phonon quantum Rabi model, approaching recent math¬ 
ematical and physical results with current quantum technolo¬ 
gies. Eurthermore, the proposed scheme provides a scalable 
quantum simulator of a complex quantum system, which is 
difficult to approach with classical numerical simulations even 
for low number of qubits, due to the large number of phonons 
involved in the dynamics. 
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Appendix A: Implementation of two-photon Dicke model with 
collective motion of N trapped ions 

In the main text, we showed how a two-photon Rabi model 
can be implemented using the vibrational degree of freedom 
of a single ion, coupled to one of its internal electronic transi¬ 
tions by means of laser-induced interactions. Here, we show 
how the N-qubit two-photon Dicke model 

n = Uja^a + ^ ^ 

n n 

can be implemented in a chain of N ions, generalizing such 
a method. The N qubits are represented by an internal elec¬ 


tronic transition of each ion, while the bosonic mode is given 
by a collective motional mode of the ion chain. The two- 
phonon interactions are induced by a bichromatic laser driv¬ 
ing with the same frequency-matching conditions used for the 
single-qubit case. The drivings can be implemented by shin¬ 
ing two longitudinal lasers coupled to the whole chain, or by 
addressing the ions individually with transversal beams. The 
former solution is much less demanding, but it may introduce 
inhomogeneities in the coupling for very large ion chains; 
the latter allows complete control over individual coupling 
strengths. 

In order to guarantee that the model of Eq. (Al) is faith¬ 
fully implemented, the bichromatic driving must not excite 
unwanted motional modes. In our proposal, the frequency of 
the red/blue drivings ujj./b satisfy the relation \oJ^/b — 0Jint\ = 
2v+5^H,, where 5^/b are small detunings that can be neglected 
for the present discussion. We recall that v is the bosonic 
mode frequency and Wmt the qubit energy spacing. To be 
definite, we take the motion of the center of mass of the ion 
chain as the relevant bosonic mode. Then, the closest collec¬ 
tive motional mode is the breathing mode [37], with frequency 
V2 = An undesired interaction between the internal 

electronic transitions and the breathing mode could appear if 
\u}r/b — 0Jint \ IS close to 1/2 or 21/2, corresponding to the first 
and second sidebands, respectively. In our case, the drivings 
are detuned by Ai = — ujint\ — ~ 0.27v from the 

first and A 2 = \u!r/b — ^int \ — ~ 1.46:^ from the sec¬ 

ond sideband. Given that the frequency 1/ is much larger than 
the coupling strength U, such detunings make those unwanted 
processes safely negligible. 

Appendix B: Properties of the wavefunctions below and above 
the collapse point 

The presence of the collapse point at g = uj/2 can be in¬ 
ferred rigorously by studying the asymptotic behavior of the 
formal solutions to the time-independent Schrodinger equa¬ 
tion Titp = Eil^- We consider now the simplest case N = 1. 
Using the representation of the model in the Bargmann space 
B of analytic functions [34], the Schrodinger equation for 
ijj{z) in the invariant subspace with generalized-parity eigen¬ 
value n = +1 reads 

#"( 2 ) + ujzijj'iz) + gz‘^'ij}{z) + ^ipiiz) = Eip{z), (Bl) 

where the prime denotes differentiation with respect to the 
complex variable z. This nonlocal linear differential equation 
of the second order, connecting the values of ip at the points 
z and iz, may be transformed to a local equation of the fourth 
order. 


ip^^\z) + [(2 — uj‘^)z‘^ + 2uj]ip"{z) -I- [4 -I- 2ujE — uj'^]zip'{z) -|- [z"* — 2a;z^ + 2 — E'^ + A‘^]ip{z) = 0, (B2) 


where we have used the abbreviations uj = ujjg, A = in the complex plane except at z = 00 , where it exhibits an 
ujqj{2g), E = E/g. Equation (B2) has no singular points 
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unramified irregular singular point of s-rank three [35]. That 
means that the so-called normal solutions have the asymptotic 
expansion 


For simplicity, we will focus on the iV = 1 case, but the pro¬ 
cedure is straightforwardly extendible to any N. We will show 
how to measure the expectation value of operators of the form 


^{z) = + ciz-i + C 2 Z-" + ...), (B3) 


exp{±m CTi (/>}crj, 


(C2) 


for z —?> oo. Functions of this type are only normalizable (and 
belong therefore to B) if the complex parameter 7 , a char¬ 
acteristic exponent of the second kind, satisfies I 7 I < 1 . In 
our case, the possible 7 ’s are the solutions of the biquadratic 
equation 


where ctij are a pair of anti-commuting Pauli matrices, 
{ai,aj} = 0, and ^ is a continuous real parameter. One 
can then reconstruct the real and imaginary parts of the 
generalized-parity operator, as a composition of observables 
in Eq. (C2) for different signs and values of i,j. 


x‘^ + x^{2-lv^) + 1 = 0. (B4) 


It follows 


^ -L 

71,2 = 2 ± 




(B5) 


For ui/2 > 1, all solutions are real. For |7i| = I74I > 1, we 
have I72I = I73I < 1. In this case, there exist normalizable 
solutions if 72 or 73 appears in Eq. (B 3 ). The condition for ab¬ 
sence of the other characteristic exponents 71 4 in the formal 
solution of Eq. (B 2 ) is the spectral condition determining the 
parameter E in the eigenvalue problem Htl’ = Eip. It follows 
that for p < a;/2, a discrete series of normalizable solutions 
to Eq. (Bl) may be found and the spectrum is therefore a pure 
point spectrum. 

On the other hand, for w /2 < 1, all 7 ^ are located on the 
unit circle with 71 = 72,73 = 74- Because, then, no normal¬ 
izable solutions of Eq. (B 2 ) exist, the spectrum of the (proba¬ 
bly self-adjoint) operator % must be continuous for p > a;/2, 
i.e. above the collapse point. The exponents 71 and 72 (73 and 
74) join at 1 (-1) for g = w/ 2 . The exponent 7=1 belongs 
to the Bargmann representation of plane waves. Indeed, the 
plane wave states (j)q{x) = ( 27 r)“^/^ expjfqa;) in the rigged 
extension of T^(K) [ 36 ], satisfying the othogonality relation 
{4‘q\4>q') = ^(9 ~ 9 O’ ^re mapped by the isomorphism I be¬ 
tween LP' (K) and B onto the functions 




-r/4 


(B 6 ) 


they correspond therefore to 7 = 1. It is yet unknown whether 
at the collapse point g = uj/2, the generalized eigenfunctions 
of % have plane wave characteristics for utq ^ 0 or which 
properties of these functions appear above this point, where 
the spectrum is unbounded from below. 


Appendix C: Generalized-parity measurement 

The generalized-parity operator, defined as If = 
(— 1 )^ cr”exp{i^n}, with n = a^a, is a non- 

Hermitian operator that can be explicitly written as the sum 
of its real and imaginary parts, 

N 

n = (-l)^(^cr"cos(|a'l’a) (Cl) 

n—1 
N 

+ sm(^a'l'a). 

n—1 


1 TT TT 

3fi(n) = --{exp(mcra;-)CT^ -I- exp(-mcra;-)cr^}, (C3) 

1 TT TT 

9(n) = -{e-xp{incF^-)ay - e-xp{-inax-)(Jy}. (C4) 

The strategy to retrieve the expectation value of observables 
in Eq. (C2) will be based on the following property of anti¬ 
commuting matrices A and B: e^Be~^ = e^^B — Be~^^. 
Based on this, the expectation value of the observables in 
Eq. (C2) can be mapped onto the expectation value of aj 
when the system has previously evolved under Hamiltonian 
H = ±nai for a time t* = (j)/2, 

{'ip\ exp{±m (Ji (j)}(Tj\'tj)) = (C5) 

where |V’(f)) = The expectation value of any 

Pauli matrix is accesible in trapped-ion setups, by fluores¬ 
cence techniques and (Jx,y by applying rotations prior to the 
measurement of The point then is how to generate the dy¬ 
namics of Hamiltonian H = ±nai. Eor that, we propose to 
implement a highly detuned simultaneous red and blue side¬ 
band interaction. 


iF = ^(o-ka1')a+e*^‘e*^-kH.c., (C 6 ) 

where (p is the phase of the laser with respect to the dipole 
moment of the ion. This Hamiltonian can be effectively ap¬ 
proximated to the second-order Hamiltonian, 

Hefi = ]{^y{2n + l)a,E^, (C7) 

when S ^ rj^Q/2. The laser phase will allow us to se¬ 
lect the sign of the Hamiltonian. Of course, one would need 
to be careful and maintain <5 in a regime where 5 v, v 
being the trapping frequency, to guarantee that higher-order 
resonances are not excited. Einally, in order to get rid of 
the undesired extra term in Hamiltonian Eq. (C7), one 
needs to implement one more evolution under the Hamilto¬ 
nian H = —(l/ 2 )Ho 77 crz. This evolution can be generated by 
means of a highly detuned carrier transition. So far, we have 
given a protocol to generate the Hamiltonian H = inuz ■ In 
order to generate Hamiltonians H = ±nay, one would need 
to modify the evolution with two local qubit rotations, 

^±inc!yt _ pia^TT/4:^±ina^t/i (C 8 ) 

Similarly, for Hamiltonian H = Enax one would need to 
perform rotations around ay. 
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